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Data collected at the Pierre Auger Observatory are used to establish an upper limit on the diffuse 
flux of tau neutrinos in the cosmic radiation. Earth-skimming v T may interact in the Earth's crust 
and produce a r lepton by means of charged-current interactions. The r lepton may emerge from the 
Earth and decay in the atmosphere to produce a nearly horizontal shower with a typical signature, a 
persistent electromagnetic component even at very large atmospheric depths. The search procedure 
to select events induced by r decays against the background of normal showers induced by cosmic 
rays is described. The method used to compute the exposure for a detector continuously growing 
with time is detailed. Systematic uncertainties in the exposure from the detector, the analysis and 
the involved physics are discussed. No r neutrino candidates have been found. For neutrinos in 
the energy range 2 x 10 17 eV < E v < 2 x 10 19 eV, assuming a diffuse spectrum of the form E„ 2 , 
data collected between 1 January 2004 and 30 April 2008 yield a 90% confidence-level upper limit 
of E'i &N Vt /AE v < 9 x 10" 8 GeV cm" 2 s" 1 sr" 1 . 



I. INTRODUCTION 

There are many efforts to search for high energy neutrinos with dedicated experiments P, S S 0> IE]- Their 
observation should open a new window to the universe since they can give information on regions that are otherwise 
hidden from observation by large amounts of matter in the field of view. Moreover, neutrinos are not deviated by 
magnetic fields and, hence, they essentially maintain the direction of their production places. The existence of ultra- 
high energy (UHE) cosmic rays of energies exceeding 10 19 eV makes it most reasonable to expect neutrino fluxes 
reaching similar energies. Although the origin of cosmic rays and their production mechanisms are still unknown 
neutrinos are expected to be produced together with the cosmic rays and also in their interactions with the background 
radiation fields during propagation 7]. Unfortunately there are still many unknowns concerning cosmic ray origin 
and neutrino fluxes remain quite uncertain. Because of their relation to cosmic ray production and transport, the 
detection of UHE neutrinos should in addition give very valuable information about cosmic ray origin. 

Models of the origin and propagation of UHECR consider the production of pions decaying into neutrinos. If 
protons or nuclei of extra-galactic origin are accelerated in extreme astrophysical environments their interaction with 
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the matter or radiation fields in the source region should yield pions which decay giving rise to neutrino fluxes. 
In addition cosmic rays interact with the background radiation when traveling over long distances resulting in a 
steepening of the spectrum around 5 x I0 19 eV. This is the Greisen-Zatsepin-Kuz'min effect consistent with the 

recently reported suppression of the cosmic ray flux above ~ 4 x 10 19 eV [HI [HJ as well as the observed anisotropy 
of the highest energy cosmic rays and a possible correlation with relatively nearby extragalactic objects [12j, |l3[. The 
GZK mechanism is a source of UHE neutrinos, in the case of protons through interactions with the Cosmic Microwave 
Background (CMB) what gives rise to the cosmological neutrinos [l4| and in the case of iron nuclei through interactions 
with infrared light that dissociates the nuclei. Alternative models, often referred to as top-down scenarios, have been 
developed although latest limits on photon fractions [15| appear to discard them as an adequate explanation of the 
UHE cosmic rays. They are based on the decay of super-massive particles into leptons and quarks. The latter 
subsequently fragment into cosmic ray protons but pions dominate the fragmentation mechanism, their decays giving 
rise to photon and neutrino fluxes. The produced neutrinos would exceed those that can be expected by the cosmic 
ray interactions with the background fields and are typically produced with harder spectra. 

Both conventional acceleration and top-down scenarios generate pions which decay to produce an electron to muon 
neutrino flavor ratio of order 1:2 while neutrinos of r flavor are heavily suppressed at production. With the discovery 
of neutrino flavor oscillations [161 ] and maximal 023 mixing, the flavor balance changes as neutrinos pro pag ate to 
Earth. After traveling cosmological distances approximately equal fluxes for each flavor are expected [13, 13- The 
idea of detecting v r induced events through the emerging r produced by neutrinos that enter the Earth just below 
the horizon, was presented for the first time in [13 . [20j . These Earth-skimming neutrinos undergo charged-current 
interactions to produce a very penetrating r lepton. When the interaction occurs sufficiently close to the Earth's 
surface the r can escape to the atmosphere and decay in flight. This would in most cases give rise to an extensive air 
shower traveling nearly horizontal and in the upward direction for an ideal spherical Earth's surface. 

The Pierre Auger Observatory (2l| has been designed to explore ultra high energy cosmic rays with unprecedented 
precision exploiting the two available techniques to detect UHE air showers, arrays of particle detectors and fluores- 
cence telescopes. It can also detect neutrinos by searching for deep inclined showers both with the surface detector [22j 
and with fluorescence telescopes [23j | . Showers resulting from t decays induced by Earth-skimming neutrinos can also 
be detected with the Pierre Auger Observatory, both with the surface and the fluorescence detectors. This channel 
has been shown to increase the possibilities for detecting neutrinos and in particular using the surface detector of the 
Pierre Auger Observatory which becomes most sensitive to neutrinos in the EeV range [24| . 

An upper limit on the diffuse flux of r neutrino of E„ &N Vir /dE v < 1.3 x 10~ 7 GeV cm~ 2 s _1 sr _1 at 90 % C.L. was 
reported in [25| using data collected between 1 January 2004 and 31 August 2007 with the surface detector of the Pierre 
Auger Observatory. The collected data were searched for r neutrino candidates applying a v T identification criterion 
that was obtained simulating Earth-skimming u T s, their interactions in the Earth, propagation of the subsequent 
t leptons and the associated showers they produce in the atmosphere. This article discusses in detail the search 
procedure to discriminate UHE Earth-skimming r neutrinos used in [25| as well as the compute of the exposure and 
the evaluation of the systematics. This article also uses an updated data sample. No candidates have been found 
in data from 1 January 2004 until 30 April 2008 and a new limit to the diffuse flux of UHE v T is presented. The 
article is organized as follows. In Section [H] the Pierre Auger Observatory is briefly described. In Section IIIII the 
needed Monte Carlo simulations are detailed. In Section [TVl the method for discriminating neutrino-induced showers 
is explained and the selection procedure is presented. In Section fV[ the computation of the exposure is reported. In 
Section [VTl the systematic uncertainties are discussed. In Section IVlIl the results from the Pierre Auger Observatory 
data for v T Earth-skimming neutrinos are shown. Finally in Section [Villi this work is summarized. 

II. THE PIERRE AUGER OBSERVATORY 

The Pierre Auger Observatory will consist of two hybrid detectors in the northern and southern hemispheres, each 
one combining an array of particle detectors and fluorescence telescopes for redundancy and calibration [26]. The 
Southern Observatory is in Malargiie, Mendoza, Argentina and its construction phase is currently completed. It covers 
3000 km 2 with regularly spaced particle detectors and with four fluorescence eyes at the perimeter that overlook the 
same area [2l[. There are plans to construct the Northern Auger Observatory in Lamar, Colorado, USA [27j • Data 
have been taken with the Southern Pierre Auger Observatory since January 2004 while it was under construction. 
The amount of data that has been accumulated for the analysis described in this article corresponds to about 1.5 
times the data that will be gathered in a whole year with the complete detector. This article will only address the 
search for Earth-skimming neutrinos with the array of particle detectors that constitutes the surface detector of the 
Southern Pierre Auger Observatory. 
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A. The Surface Array of the Pierre Auger Observatory 

The Southern surface detector (SD) consists of 1600 Cherenkov water tanks (3.6 m diameter and 1.2 m high) 
arranged in a triangular grid with 1.5 km spacing between them, covering an almost flat surface, at an approximate 
altitude of 1400 m above sea level. Each tank is a polyethylene tank internally coated with a diffusive Tyvek liner 
filled with 12 tons of purified water. The top surface has three Photo Multiplier Tubes (PMTs) in optical contact 
with the water in the tank. The PMT signals are sampled by 40 MHz Flash Analog Digital Converters (FADC). Each 
tank is regularly monitored and calibrated in units of Vertical Equivalent Muons (VEM) corresponding to the signal 
produced by a /i traversing the tank vertically [2Sj ] - The system transmits information by conventional radio links 
to the Central Data Acquisition System (CDAS) located in Malargiie. The PMTs, a local processor, a GPS receiver 
and the radio system are powered by batteries with solar panels. Once installed, the local stations work continuously 
without external intervention. 

The local trigger at the level of an individual Cherenkov tank (second order or T2 trigger) is the logical OR of two 
conditions: either a given threshold signal (1.75 VEM) is passed in at least one time bin of the FADC trace, or a 
somewhat lower threshold (0.2 VEM) is passed at least in 13 bins within a 3 fxs time window (120 bins) [29(. The 
latter condition, the so-called Time over Threshold (ToT), is designed to select broad signals in time, characteristic 
of the early stages of the development of an extensive air shower (EAS). The data acquisition system receives the 
local triggers and builds a global trigger requesting a relatively compact configuration of 3 local stations compatible 
in time, each satisfying the ToT trigger, or 4 triggered stations with any T2 trigger (a third level or T3 trigger) [3(| • 
With the complete array, the global T3 trigger rate will be about 3 events per minute, one third being actual shower 
events at energies above 3 x 10 17 eV. 

B. The Data Sample 

The SD has been taking data in a stable manner since January 2004 [3l[ . Meanwhile the array has been growing 
and the number of deployed stations has increased from 120 to 1600 during the period analyzed in this article. 

The analysis reported here is restricted to selected periods in order to eliminate inevitable problems associated to the 
construction phase, typically in the data acquisition, in the communication system and due to hardware instabilities. 
To ensure the quality of the data, we have analyzed the arrival time of the events under the reasonable hypothesis that 
the rate of physics events recorded by the detector (after proper size normalization) is independent of time. Given the 
large aperture of the SD and the level at which anisotropies could exist on the sky [32 | this approximation is, from 
this point of view, well justified. Assuming a constant rate A of physics events, the probability P of the time interval 
t between two consecutive events to be larger than T is given by: 

P(t > T) = e- XT (1) 

where the value of A is the mean rate of the recorded events normalized to the detector size. 

Consecutive events for which P is below a certain threshold value P cu t are assumed to belong to periods with 
problems in the data acquisition and are used to define the bad periods to be rejected. The procedure will reject a 
good period with probability P C ut together with the eventual bad ones. So in principle one would like P cu t to be as 
small as possible. The choice of P cu t is made by finding where the distribution of probability of the time interval 
between two consecutive events differs from being flat. The numerical value of P cu t was found to be ~10 -5 , which 
allows us to reject only a small fraction of good periods while removing the periods that lead to a non flat probability. 
When data taking began, the bad periods were of the order of 10% of the operating time but by the end of the time 
period considered in this work we were typically below 1% 33]. 

Once the bad periods have been removed, the events that have passed the third level trigger [3(| from January 2004 
until April 2008 constitute the data sample used in this paper. 

III. END TO END SIMULATION CHAIN 

In order to obtain a flux or a flux limit from the data, the detector neutrino showers have to be searched with 
a selection criterion and the exposure of the detector must be accordingly computed. Both the criteria to identify 
neutrino induced showers and the computation of the exposure to v T are based on Monte Carlo techniques. Three 
separate simulations can be identified. Firstly a dedicated simulation that deals with the neutrinos entering the Earth 
and the r leptons that exit. A second simulation involves the r decay in flight and the development of an up- going 
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atmospheric shower. Finally a simulation of the tank response to the through-going particles is performed to convert 
the particles at ground level obtained in the shower simulation to an actual detector signal. 

A. Earth-Skimming neutrinos 

As Earth-skimming v T s penetrate the Earth they interact to produce r leptons that suffer energy loss but can 
escape the Earth and decay in the atmosphere. As a result the incoming neutrinos give rise to an emerging r flux 
which depends on the depth of matter traversed (for a spherical Earth as assumed here, it depends only on the nadir 
angle). The decay of the r lepton in the atmosphere induces an extensive air shower (EAS) that can trigger the 
SD. The efficiency of this conversion plays a key role in the calculation of the detector exposure. The r flux has 
been computed using simulation techniques that take into account the coupled interplay between the r and the v T 
fluxes as they traverse large matter depths through charged current (CC) weak interactions and through r decay. 
Energy losses induced by neutral current (NC) interactions for both particles are taken into account as a stochastic 
process. The energy losses through bremsstrahlung, pair production and nuclear interactions for the r lepton are 
applied continuously, which, at the level of accuracy we need for this work, is a reasonably good approximation [34j . 

Propagation of particles through matter is performed in small depth steps. At each step the particles are followed 
and the probability for interaction and decay (in the case of the r) are evaluated taking into account the particle 
energy. The chain starts with an incident v T which may interact by CC or NC. When the former occurs, a r lepton is 
generated, and its energy is selected taking into account the ^-distribution of the interaction, where y is the fraction 
of the v T energy transferred to the nucleon in the laboratory frame. If the v r interacts through a NC its energy 
is computed taking the y-distribution into account. Once a r lepton is produced, it can undergo energy loss, weak 
interactions both neutral and charged and decay. In the case of CC interaction or decay a new v T is produced which 
regenerates the v r flux that is propagated further. Finally, if a r lepton emerges from Earth, its energy, direction and 
decay position are stored and used as an input for the simulation of atmospheric showers induced by r leptons. 

For the relevant depths inside the Earth where r leptons can be produced and reach the surface (few km) a 
homogeneous density of 2.65 g cm -3 can be assumed. Parameterisations of the cross section for weak interactions 
and for the ^-distributions at very high energy are used. The cross section for CC interactions is taken fro m 13511 and 
the y-distribution from j36j. For NC interactions, the cross section is assumed to be 0.4 that of the CC [37| . The 
energy losses for t leptons are parameterised following case III in [38( , which gives the best representation of Monte 
Carlo simulation. 



B. Extensive Air Showers in the atmosphere 

The t decays in the atmosphere give rise to secondaries that may initiate an EAS that can trigger the SD. The 
decay mode has been simulated using the TAUOLA package version 2.4 l39ll to obtain the type of the secondaries and 
their energies which are subsequently injected in Aires (version 2.6.0) |40l] with SIBYLL 2.1 [41] as a model for the 
hadronic interactions at high energy. Showers induced by up-going rs with energies from log(E T /eV) = 17 to 20.5 
in steps of 0.5 have been simulated at zenith angles ranging between 90.1° and 95.9° in steps of 0.01 rad and at an 
altitude above the Pierre Auger Observatory that ranges from m to 2500 m in 100 m steps. Ten showers have been 
generated for each combination of energy, zenith angle and altitude, which leads to a total of 20 000 showers. 

The extremely large amount of particles involved in an ~ EeV shower makes it impractical to follow all the 
secondaries. The current simulation packages include a statistical sampling algorithm based on the thinning algorithm 
originally introduced in (42J. Only a small representative fraction of the total number of particles is propagated. 
Statistical weights are assigned to sampled particles in order to compensate for the rejected ones. 

C. Detector response 

The first step in the detector simulation is to obtain the particles reaching each tank from the sampled particles 
produced in the simulation of the EAS. A re-sampling algorithm is necessary to convert the output of the program 
to the expected number of particles that enter a SD station. This is done averaging over an area around the station 
that is large enough to avoid unphysical fluctuations from the thinning procedure, and at the same time small enough 
to avoid large differences in the density and average properties of particles in different places on the area [43[ . Each 
particle reaching the station is injected inside the tank, and a detailed simulation is performed to obtain the light 
hitting the PMTs as a function of time. The simulated FADC traces are obtained as the superposition of the signal of 
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all individual particles entering the tank accounting for their arrival time. Finally both the local and central trigger 
algorithms are applied and the event is stored in the same format as data f44| . 

At the highest simulated values of incident angles or altitudes where the r decays, none of the simulated showers at 
any of the simulated energies fulfills the central trigger conditions. This is taken as a clear indication that a complete 
sample of showers has been produced without introducing any bias and that it therefore correctly represents the 
characteristic r showers that could trigger the SD detector. 

IV. DISCRIMINATION OF NEUTRINO-INDUCED SHOWERS 

A. Neutrino signature: inclined showers in the early stages 

UHE particles interacting in the atmosphere give rise to a shower with an electromagnetic component reaching its 
maximal development after a depth of the order of 800 g cm -2 and extinguishing gradually within the next 1000 
g cm~ 2 . After roughly a couple of vertical atmospheric depths only high energy muons survive. In the first stages 
of development, while the electromagnetic component develops, the time spread of the particles in the shower front 
is large (~ (is). When the shower becomes old, most of the particles in the shower front, the high energy muons, 
arrive in a short time window ( ~ 100 ns ). As a consequence very inclined showers induced by protons or nuclei 
(or possibly photons) in the upper atmosphere reach the ground as a thin and flat front of muons accompanied by 
an electromagnetic halo, which is produced by bremsstrahlung, pair production and muon decays, and has a time 
structure very similar to that of the muons. On the other hand, if a shower is induced by a particle that interacts 
deep in the atmosphere (a deep neutrino interaction in air, or a tau decay), its electromagnetic component could hit 
the ground and give a distinct broad signal in time. The signal in each station of the SD is digitized using FADCs, 
allowing us to unambiguously distinguish the narrow signals from the broad ones and thus to discriminate stations 
hit by an EAS in the early stages of development or by an old EAS. This is illustrated in figure [1] where we show 
FADC traces from two different real events. The FADC trace taken from the shower with a zenith angle of 22° is 
representative of an EAS in the early stages while the other is representative of an old EAS. 

B. Identification of neutrino candidates 

The identification of showers induced by Earth-skimming r neutrinos implies searching for very inclined (quasi- 
horizontal) showers in an early stage of development. Broad signals, which are characteristic as long as the electro- 
magnetic component still develops, produce a ToT local trigger (see section Hi A [I . A bunch of muons from cosmic ray 
showers can produce high amplitude signals extended in time or two independent muons can arrive inside the given 
time interval. Both would also produce a ToT local trigger which is not associated to the presence of electromagnetic 
component from a neutrino shower. To get rid of them, a further requirement is made to the signals in order to 
filter out these backgrounds. First a cleaning of the FADC trace is done to remove segments of the trace that could 
be generated by an accidental muon arriving closely before or after the shower front. Segments of the FADC trace 
are defined by neighbour bins above 0.2 VEM, allowing gaps of up to 20 bins and only the segment with the largest 
signal is kept. An offline ToT is defined by requiring that the signal after cleaning (the segment with largest signal) 
of the FADC trace has at least 13 bins above the low threshold (0.2 VEM) and the ratio of the integrated signal over 
the peak height exceeds by a factor 1.4 the average ratio observed in signals of isolated particles (as defined in the 
calibration procedure [HI]). The central trigger conditions are applied only to stations that fulfill the offline ToT. Still 
a small number of nucleonic showers with a large number of triggered tanks may have a subsample of stations that 
satisfy this condition even if in all the other stations the signal is not broad at all. In order to reject such events, 
at least 60% of the triggered stations are required to fulfill the offline ToT. After this selection procedure, an almost 
pure sample of showers reaching the ground at their early stages is selected. 

Once the criterion for young showers is established a second criterion must be used to select very inclined showers 
as expected from Earth-skimming neutrino interactions. The devised method uses two variables associated to the 
footprint that the triggered tanks of the event leave on the ground and the apparent speed with which the signal 
moves across the array. Firstly a symmetric tensor is built using the station signals included in the central trigger 
and their ground positions (analogous to the tensor of inertia of a mass distribution, see Eq. J2])). 

S = 'Ei-Si {X} = T,iSiXi/S (Y) = ^SiUi/S 

I xx = X iSi (xi - {X)f/S Iyy = ViSiivi - {Y)f/S 

I xy = I yx = Xfstixi - (X))( yi - {¥)) (2) 
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Energy of shower ~ 5 EeV 
Distance to shower axis ~ 1 .0 km 
Zenith angle ~ 22° (early stage) 
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FIG. 1: FADC traces of stations at 1 km from the shower core for two real showers of 5 EeV. Top panel: early stages of 
development (8 ~ 22°); bottom: old extensive air shower (0 ~ 80°). 



where Si is the signal in VEM for each station; (xi, j/j) are the coordinates of each station; and is the sum over the 
stations. 

The corresponding major and minor axes are used to define a characteristic length and a width of the pattern as 
the square root of the eigenvalues of the symmetric tensor (see Eq. Q). 



Ixx + lyy + J (Ixx lyy) 2 + ^xy 

length = — 



Ixx + lyy \ {Ixx lyy) 2 + 

width 2 = ¥_ (3) 

Secondly for each pair of tanks a ground speed is defined as dij/\Atij\, where dij is the distance between 

the tanks (projected onto the major axis) and |Atjj-| is the difference between the start times of their signals (figure 
[2]). Quasi-horizontal showers have an elongated shape (characterized by a large value of length/ width) and they have 
ground speeds tightly concentrated around the speed of light c. 

In figure[3]the distributions of the two discriminating variables are shown for real events and simulated tau showers. 
Based on the comparison between MC simulations and data collected during November and December 2004, which is 
less than 1% of the used data sample, the following cuts were fixed to select Earth-skimming tau neutrino candidates: 




FIG. 2: Schematic view of the footprint of a shower on the SD array. Each circle represents the position of a station, and their 
sizes are proportional to the station signal. 




FIG. 3: Distribution of variables used to discriminate very inclined showers for an incident E~ 2 v T flux (histogram), and for 
real events collected during November and December 2004 passing the early stage (see text) selection (points). Left panel: 
length/width; middle: average of the ground speed between pairs of stations; right: r.m.s. of the ground speeds. 

• length / width > 5 

• 0.29 m ns _1 < average ground speed < 0.31 m ns _1 

• r.m.s. (ground speed) < 0.08 m ns _1 

where the average ground speed and its dispersion are computed using only stations for which \di t j \ is larger than 1000 
m. 

Since the selection criteria relies on the shower footprint, we need to guarantee that a representative fraction of 
the event is detected with the SD. For this purpose the closest station to the center of the footprint (values (X) and 
(Y) defined in Eq. ^ is required to be surrounded by at least five working stations at the time of occurrence of the 
event. Hence, events at the edges of the array with small detected fraction of the footprint typically do not fulfill the 
selection criteria. This procedure is simple and robust. It can be applied to any footprint and does not require any 
global reconstruction. 

V. NEUTRINO EXPOSURE OF THE SURFACE DETECTOR OF THE PIERRE AUGER 

OBSERVATORY 

The next step in the calculation is to compute the exposure of the SD of the Pierre Auger Observatory to showers 
induced by UHE v T . Each simulated Earth-skimming v T event has to be tracked from the injection up to its identi- 
fication through the defined selection cuts. The number of identified events is computed from the simulations of the 
EAS initiated by the secondaries in the tau decay, and from the detector response to them. For a fixed energy of the 
r (E T ), there is effectively only one relevant parameter determining the efficiency of trigger and identification: the 
altitude of the shower center (h c ). This is conveniently defined as the altitude of the shower axis at a distance of 10 
km away from the r decay point along the shower axis (see figure 2]). For the shower energies relevant in this analysis, 
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h c is very close to the altitude at which a horizontal shower has the largest lateral extension and is thus capable of 
producing the largest footprint at ground [24| . 




FIG. 4: Geometry of the induced r shower with the definition of the parameters, h c and O, involved on the exposure computation 
(see text). Angles and distances are not proportional. They have been exaggerated to help the readability of the figure. 

For a SD that covers a surface A, the aperture for a given neutrino energy (E v ) can be expressed as follow: 

Ap(E v ) = [ AcosQdVL f dE T [ dh c ( 

J Jo Jo V dE T dh c 

fir/2 pE v poo ( cfin (F )\ 

= 2tt / AcosOsinOdO / dE T / dh c , Uff 

J*/2+ am JO JO V dE T dh c J 

= ,«,V(kf*(«) (/ , (4) 

where d 2 p T /(dE T dh c ) is the differential probability of an emerging r as a function of energy and altitude for a hxed 
incident v T energy, that can be easily obtained folding the simulations described in section IIII Al for the emerging rs 
with the tau decay probability as a function of flight distance. £// is the probability to identify a r (including the 
trigger efficiency), that is assumed to depend only on E T and h c . The integral in is done from 7r/2+a m rad (a m =0.1 
rad) to 7r/2 rad since an incident v T with a greater angle has no chance to produce an emerging r that produces an 
observable shower at ground level. The latter integration can be performed by the Monte Carlo technique as described 
in section Hill leading to: 




Ap{E v ) = ttW a m ^ eff ^ " ^ ] (5) 

where N S i m is the number of simulated events. In figure [5l the trigger and identification efficiencies for an ideal (no 
holes and no malfunctioning stations) and infinite array are shown as given by the simulated EAS described in section 
IIII Bl They have been calculated by throwing once each simulated EAS on the detector array with a random core 
position. The maximum efficiency that can be reached is 82.6 % due to the /i channel decay [451 ]. This decay mode 
does not produce a detectable shower neglecting the possibility of hard muon bremsstrahlung or pair production near 
the detector which should have a negligible effect on the final limit. The identification efficiency depends smoothly 
on E T and h c , and hence it can be safely interpolated. 
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During the period of data taking considered in this work, the SD of the Pierre Auger Observatory has been growing 
continuously. It is of course mandatory to take into account this evolution, as well as the instabilities of each station. 
Therefore Eq. ((5]) is not valid to compute the actual SD array. Instead the following expression can be used: 



Ap(E v ,t) = 



(6) 



-it sin 2 a, 



dE T I dh c 



d 2 p T 
dE T dh l 



dxdy eff(E T ,h Q x,y,Aconf(t)) 



where e// now also depends on the position of the shower in the array (x,y), and on the instantaneous configuration 
of the array at time t denoted here as Aconf(t)- The integral over the area S includes the whole SD array. 
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FIG. 5: Trigger (open dots) and identification (closed dots) efficiency as a function of the height above ground of the shower at 
10 km from the decay point. The efficiency is shown for MC showers induced by ts with energy of 0.1 (top-left), 1 (top-right), 
10 (bottom-left) and 100 (bottom-right) EeV merging all zenith angles. 

Hence, the total exposure during the considered period of data taking is the time integration of the instantaneous 
aperture given by: 



Exp = I dtAp{E v ,t) =7rsin 2 



dE T 



dh c ( d2pT B T 
u c \ dE T dh c 



B T (E T ,h c ) = dt / dxdyeff(E T ,h c ,x,y,A Co nf(t)) 



(7) 



The exposure is computed by Monte Carlo in two independent steps. First, the integrals in t and (x,y) are 
computed using the simulations of the EAS and the detector. The number of working stations and their status are 
monitored every second allowing us to know with very good accuracy the instantaneous SD configuration [46]. For 
each simulated EAS, several random times from January 2004 until April 2008 excluding the rejected periods are 
selected. The number of random times is selected in a monthly base to ensure a statistical precision on the exposure 
at 1% level. For each time, the evaluation of the identification efficiency is done for any position of the shower in the 
SD array. The average over all showers with the same E T and h c gives the integral in time and area of e//, allowing 
one to compute B T (E T , h c ). The second step computes the integral in h c and E T as in the case of a perfect array. 
The estimated uncertainty of this method given the Monte Carlo simulations is below 3%. The accumulated exposure 
is shown in figure [SJ It corresponds to an equivalent time of about 1.5 years of the complete SD array (1600 water 
tanks). 
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FIG. 6: Earth- Skimming neutrino exposure of the Pierre Auger Observatory accumulated from January 2004 until April 2008. 
All the identification cuts described in setion lTV Bl are taken into account. 

VI. SYSTEMATIC UNCERTAINTIES 

Several sources of systematic uncertainty have been carefully considered. They are addressed below. We have 
chosen as a reference the aperture calculated with v cross section from Ref . (35j , the parameterization of the energy 
losses from Ref. [38| , an uniform random distribution for the r polarization, a spherical model of the Earth and the 
SIBYLL [4l[ hadronic model in combination with AIRES shower simulator [4(| ■ The systematic uncertainties in this 
section are all quoted with respect to this aperture and therefore in general asymmetric. Moreover, to be able to 
quote a range for the systematic uncertainties independently of the energy, an E~ 2 incident flux of neutrinos has been 
assumed. 

Firstly, the location of the Pierre Auger Observatory is close to the Andes and not very far away from the Pacific. The 
actual topography of the Pierre Auger Observatory can be taken into account by a detailed Monte Carlo simulation [t^ ]. 
The effect of the Andes on the expected event rate has been studied with the aid of a digital elevation map available 
from the Consortium for Spatial Information [48]. The number of detected events decreases by 18 % if Andes are 
neglected. 

There is quite some level of uncertainty in EAS simulation because accelerator data have to be extrapolated to the 
shower energies under discussion. However these uncertainties are not expected to have a large effect on the final 
result since the electromagnetic component of the shower, which is the most relevant part for neutrino identification, 
is believed to be better reproduced by simulations. Shower simulations have been done with two hadronic models 
(QGSJET [49[ and SIBYLL [4l[) and passed through two different detector simulations. Based on that, systematic 
uncertainties of are quoted as due to the Monte Carlo simulation of both the EAS and the detector, the former 

being the main contribution 1 . The simulations of the interactions inside the Earth have been extensively checked 
by comparison with an analytical calculation [5lT ] , an iterative solution of the transport equations [52| , and several 
independent simulations. The uncertainty associated to the simulation process itself is expected to be below the 5 % 
level. 

Monte Carlo simulations also make use of several physical magnitudes that have not been experimentally measured 



1 Currently, QGSJET and SIBYLL are the only hadronic models available to be used with the used EAS simulation package. Other 
recently introduced models like EPOS [50| are not yet available to test their effect. Although, in the case of EPOS, due to the large 
number of muons, and the flatter lateral distribution the trigger efficiency should be larger and in this respect our limit should be 
conservative. 
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at the relevant energy range, namely: the v cross-section, the r energy losses and the polarisation of the r. All of them 
can be computed in the framework of the Standard Model of particle physics using the parton distribution functions 
(PDFs). 1 _ 

The allowed range for the v cross section due to uncertainties in the PDFs has been studied in [35j and includes 
both the effects of the experimental uncertainties on the PDFs fitted to ZEUS and fixed target data evolved at next- 
to-leading-order [53|], as well as theoretical uncertainties in the implementation of heavy quark masses on the PDF 
evolution. For the purpose of this study, the ZEUS PDFs and their uncertainties were recalculated (3f| using the 
DGLAP equations throughout the relevant kinematic range (down to x ~ 10~ 12 ). This leads to a systematic 
uncertainty for the number of v expected to be detected by the SD of the Pierre Auger Observatory. 

The decay of the r lepton plays a key role on the whole Monte Carlo simulation. Both the branching ratios of 
the different decay modes and the energy distribution among the products are important. The latter depends on 
the r polarisation, which in turn depends on the PDFs. The most and least favorable cases in the range of possible 
polarisations (helicity ±1) have been used to estimate the uncertainty associated to it. The use of the extreme cases 
of polarisation of the r will not produce more than t\o% differences on the exposure. 

Finally, energy losses include r bremsstrahlung (BS) and pair production (PP) as well as nuclear interactions. The 
contributions from BS and PP can be accurately rescaled from the values for muons [45ll54T|. The nuclear contribution 
comes from the photo-nuclear cross-section and it is much more uncertain. The differential photo-nuclear cross-section 
as a function of the PDFs has been given in [EH, [56[ . There exist estimates of the tau energy losses for the relevant 
energy range based on them [23|, [H, HH, HH, [57J • Different calculations of the energy losses may lead up to 
systematic uncertainties in the exposure. 

In figure [7] the exposure for 1 year of the SD array with 1600 water tanks is shown in the most and least favorable 
cases of the systematic uncertainties previously discussed. The systematic uncertainties do not have the same effect 
for all v T energies. The importance of each different contribution to the global systematic uncertainty in the exposure 
is neither the same at all v T energies. At low energies (~1 EeV) the r polarization, the v cross section and the r energy 
losses dominate. At higher energies those contributions become smaller and other increase. The latter comes from 
neglecting the mountains, the effect of which increases with the v T energy. The small depth traversed becomes more 
relevant due to the larger cross section. The former is due to the contribution from the v cross section uncertainties. 
A larger cross section increases the interaction probability for the neutrinos but reduces the solid angle due to the 
flux absorption in the Earth. This makes the uncertainty in the cross section to contribute mainly around 1 EeV. 

The effect of the systematic uncertainties on the expected rate of identified v T will depend on the shape of the 
actual incident v T flux. The effect is almost the same either for GZK-like fluxes or for E~ 2 fluxes, giving a factor 
~ 3 for the systematic uncertainty in either case (see Table H}. Moreover, the energy dependent effect also produces 
differences on the energy range where most of the identified v T are expected. For instance, if an E~ 2 flux is assumed 
the energy range where 90% of the events are expected changes from 0.22-23 EeV in the least favorable scenario to 
0.20-26 EeV in the most favorable one. 



Source 


factor 


EAS Simulations 


1.30 


Topography 


1.18 


Tau Polarisation 


1.30 


Cross Section 


1.15 


Energy losses 


1.40 


Total 


2.9 



TABLE I: Ratio of expected number of v T for either GZK-like or for E 2 incident spectra in the most and least favorable 
scenarios for each source of systematic uncertainties. 



The relevant range of PDFs involved in both the v T and the r photo-nuclear cross-sections includes combinations 
of Bjorken-a; and Q 2 where no experimental data is available. Only extrapolations that follow the behavior observed 
in the regions with experimental data have been considered. Different extrapolations to low x and high Q 2 would 
lead to a wide range of values for the v cross-section as well as for the r energy losses. The systematic uncertainties 
due to this have not been included in the quoted systematics. Possible large v cross-sections have not been taken into 
account either. 
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FIG. 7: The Pierre Auger Observatory exposure for one year of a full SD array data. Solid lines bracket the total systematic 
uncertainty, while the dotted ones only bracket the allowed range of polarisations, v cross-sections and energy loss uncertainties. 



VII. SEARCH RESULTS AND NEUTRINO LIMIT 



The measurement of the spectral shape or energy dependent upper limits on v T are out of the reach of Pierre 
Auger Observatory scope for mainly three reasons. Firstly the generic energy reconstruction algorithms developed for 
conventional nucleonic showers do not work if the position of the shower axis is not determined (here, a knowledge of 
the altitude at which the shower is produced would be needed); secondly, the fraction of r energy contributing to the 
EAS depends on the decay mode; and finally, the energy transferred from the incident v r to the emerging r is not 
known. At best an approximate lower bound of the initial v T energy could be obtained. 

The data have been searched for neutrino candidates over the analysed data period and there are no events fulfilling 
the selection cuts. The events failing to pass only one of them have a quite different distribution for the discriminant 
parameter than the one of the simulated r showers (see figure [8]) . In Table HH the number of events surviving after 
successive cuts for real data as well as the efficiencies for simulated v T s are shown. The huge reduction of events after 
selecting very inclined showers (Elongated footprint and Ground speed) reaching the ground in early stages (Young 
showers) is expected for showers induced by protons or nuclei (see Section TlV Ap . The expected number of events 
surviving the combination of those cuts due to detector effects is also compatible with the observed discriminating 
power. Based on that, a limit for an injected spectrum K ■ Q(E) with a known shape Q(E) can be derived. The 90% 
confidence-level (CL), for candidates and no background expected [58], on the value of K is: 



K9 ° ~ f*(E).Exp(E)dE (8) 

where Exp is interpolated from Table HTT1 

For an injected diffuse flux of v T dN/dE = K ■ E~ 2 , the 90% CL limit is K go = 6±^ x 10~ 8 GeV cur 2 s" 1 sr~\ 
where the uncertainties come from the systematics discussed in section IVI1 The bound is obtained for the energy 
range 2 x 10 17 - 2 x 10 19 eV, with a systematic uncertainty of about 15%, over which 90% of the events can be expected 
for an E~ 2 flux. 

In figure [9] the limit for the most pessimistic scenario of systematic uncertainties is shown. It improves by a factor 
~ 3 in the most optimistic scenario (dotted line). Flux limits given by other experiments are also shown (divided by 
3 if they are limits to all flavours to be able to compare): AMANDA [H H?|, Baikal [gj, RICE [g^ (rescaled at 
90 % CL), HiRes [H Hi], ANITA-lite [§3, ANITA [l36[, GLUE [g7| and FORTE For some of them the limits 
are given for an E~ 2 flux (integrated format), while for others the flux limit is given as 2.3/Exposure-i?^ (differential 
format). The limit from the ANITA-lite experiment and the Pierre Auger Observatory are given in both formats for 
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FIG. 8: Distribution of discriminating variables for an incident E~ 2 v T flux (histogram), and for the real events (points). 
Events that survive all the selection cuts except the one for the shown variable are used. The vertical lines indicate the values 
to discriminate. Left panel: fraction of stations with a ToT signal; middle: length/width ; right: average of the ground speed 
between pairs of stations. 



Selection Requirement 


MC Efficiency 


Number of real events 


Initial sample 


1.00 


3.97xlO b 


Young showers 


0.88 


6.68xl0 5 


Elongated footprint 


0.87 


8.37 xlO 3 


Ground speed ~ c 


0.84 





Contained footprint 


0.76 






TABLE II: Number of events passing the successive selection cuts. The Monte Carlo efficiency corresponds to identification 
efficiencies for ^s of energy 1 EeV that trigger the SD detector. Data have been collected from January 2004 until April 2008. 



comparison. The differential format demonstrates explicitly that the sensitivity of the Pierre Auger Observatory to 
Earth-skimming neutrinos peaks in a narrow energy range close to where the GZK neutrinos are expected. 

The energy range of the v T s explored with the Auger Observatory with this channel is very well suited to search 
for the diffuse flux of ^s that are produced by the GZK effect. 



VIII. SUMMARY AND PROSPECTS 



The dataset collected during the construction phase of the surface detector of the Pierre Auger Observatory from 
January 2004 until April 2008, is used to present an upper limit to the diffuse flux of v T . The Earth-skimming 
technique together with the configuration of the surface detector gives the best sensitivity currently available around 
a few EeV, which is the most relevant energy to explore the predicted fluxes of GZK neutrinos. However in the worst 
case of systematic uncertainties, the limit presented here is still higher by about one order of magnitude than GZK 
neutrino predictions. The Pierre Auger Observatory will keep taking data for about 20 years over which the bound 
will improve by more than an order of magnitude if no neutrino candidate is found. 
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xlO 16 
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xlO 16 
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xlO 17 
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xlO 16 


20.5 
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xlO 17 
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xlO 16 



TABLE III: Exposure of the SD of the Pierre Auger Observatory from January 2004 until April 2008. In the columns labeled 
highest and lowest we give the values of the exposure for the most optimistic and most pessimistic cases of the systematic 
uncertainties. 




10 14 10 16 10 18 10 20 10 22 10 24 10 26 

Neutrino Energy [eV] 

FIG. 9: Limits at 90 % CL for each flavor of diffuse UHE neutrino fluxes assuming a proportion of flavors of 1:1:1 due to 
neutrino oscillations. The Auger limits are given using the most pessimistic case of the systematics (solid lines). For the 
integrated format, the limit that would be obtained in the most optimistic scenario of systematics is also shown (dashed line). 
See text for the references to the other experimental limits. The shaded area corresponds to the allowed region of expected 
GZK neutrino fluxes computed under different assumptions [69L Frol. FtD. l72| . although predictions almost 1 order of magnitude 
lower and higher exist. 
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